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We report a novel algorithm of constructing linear and nonlinear potentials in the two-dimensional 
Gross-Pitaevskii equation subject to given boundary conditions, which allow for exact analytic so- 
lutions. The obtained solutions represent superfluid flows in inhomogeneous Bose-Einstein conden- 
sates. The method is based on the combination of the similarity reduction of the two-dimensional 
Gross-Pitaevskii equation to the one-dimensional nonlinear Schrodinger equation, the latter allowing 
for exact solutions, with the conformal mapping of the given domain, where the flow is considered, 
to a half-space. The stability of the obtained flows is addressed. A number of stable and physically 
relevant examples are described. 
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I. INTRODUCTION 

Nowadays one observes rapidly increasing interest in 
studying nonhnear Schrodinger (NLS) equations with in- 
homogeneous coefficients and, in particular, in obtain- 
ing their exact analytical solutions for the physically 
relevant statements. Starting with the first results on 
integrable inhomogeneous models p^, this activity re- 
ceived further development due to its relevance in non- 
linear optics [2] and in the mean-field theory of Bose- 
Einstein condensates (BECs) d, Q • A possibility of con- 
structing exact solutions was also reported for NLS equa- 
tions with inhomogeneous complex- valued coefficients [5| . 
However, there are two main limitations of the presently 
available results. First, most of them were obtained 
for one-dimensional or quasi-one-dimensional statements 
and only a few results on multidimensional problems were 
reported, so far 0, Q. Second, the most of the models 
allowing for construction of exact solutions were posed in 
the infinite domains. The first suggestion of an algorithm 
for constructing exact solutions of the one-dimensional 
(ID) NLS equation on a half-line, modeling a BEC inter- 
acting with a rigid surface, was recently reported in 0. 

Here we show that the requirements of one- 
dimensionality and unboundness of the domain can be 
removed and exact analytical solutions can be obtained 
for models defined on bounded 2D domains and de- 
scribed by the NLS equation with inhomogeneous linear, 
T4xt(i*), and nonlinear, ^(r), potentials (the both being 
real- valued functions of the spatial coordinates). More- 
over, some of the reported solutions are found to be sta- 
ble, and thus having particular physical relevance. 

The paper is organized as follows. In Sec. II, we 
present the 2D physical model and give the conformal 



mapping to reduce the 2D physical model subject to 
given boundary conditions to nonlinear ordinary differ- 
ential equation solved. In Sec. HI, we concentrate on 
the two simplest representative examples to illustrate our 
novel method. Sec. IV is devoted to the numerical sim- 
ulations for the solutions obtained in the domain D2. In 
Sec. V, we study the generalization of the conformal 
mapping and present exact solutions of the 2D physi- 
cal model. Finally, the outcomes are summarized in the 
Conclusion. 



II. THE PHYSICAL MODEL AND 
CONFORMAL MAPPING 

To be specific, we deal with the 2D nonlinear physical 
model 



■^V2 + T4xt(r)+5(r)|*(r,t)|^ 



*(r,t),(l) 
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where r = {x^y) G C M^, I) is an open domain, and 
V = (dx^dy). We are particularly interested in applica- 
tions of our results to BEC flows, where ^(r, t) is the 
macroscopic wavefunction and modelJT]) is also termed, 
the Gross-Pitaevskii (GP) equation [8[. We explore the 
flexibility of potentials in the BEC applications, i.e., pos- 
sibility of manipulating them by external electric and/or 
magnetic fields for the sake of creation of desirable spatial 
configurations for the linear potential and for the scatter- 
ing length of the two-body interactions (the latter per- 
formed through the Feshbach resonance technique [9[). 
We also notice that the model ([1]) has also direct rele- 
vance to the mean-field theory of exciton-polariton con- 
densates [lo|- There on the one hand, the 2D statement, 
i.e. the statement considered in this paper, is the most 
typical one. On the other hand applying the external 
pump from the free edges of a specimen one can create 
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different kinds of nonzero conditions (nonzero currents, 
as required bellow in the present paper). 

We concentrate on a BEC in a domain D bounded by 
impenetrable walls. Respectively Eq. ^ will be supplied 
by the zero conditions given at the boundary of the do- 
main which we denote as dD^ i.e., we impose ^(r) = 
for ah r G dD. 

Our goal is to find an algorithm allowing for system- 
atic constructions of the linear, lext(r), and nonlinear, 
^(r), potentials, for which the formulated Dirichlet prob- 
lem allows for exact analytical solutions. To this aim we 
assume that there exists a complex analytic function 

C(r) = r,(r) + i^{v) = f{z) (2) 

of the complex variable z = x-\-iy G C, which provides the 
conformal mapping of the contour dD to the imaginary 
axis of C{t)^ i.e., to 77(r) = 0, such that the domain D 
is mapped into the right half-plane in terms of the new 
variables {r]^(f): r]{r) > 0. Due to analyticity of the 
mapping C(r) the Cauchy-Riemann equations on the pair 
of real- valued functions 7^(r) and (p{r) hold: 

da^rjir) = dycpir), dyr](v) = -d^^ifiv). (3) 

Respectively, the following constraints 

V^7?(r) = 0, VV(r) = 0, V r]{Y) -V Lp{v) = (4) 

are verified, as well. In other words, 77(r) and (p{r) are 
both 2D harmonic functions with orthogonal gradients. 
Moreover, on the basis of the Cauchy-Riemann equations 
we find the relation |V7^(r)p = |V(/?(r)p. 

We restrict the consideration to the linear and nonlin- 
ear potentials which can be represented in terms of 77 (r) 
as follows 

yext(r)^-||V7?(r)|2, 5(r)^||Vr,(r)|2 (5) 

where e and Q are real parameters. Without loss of 
generality, we choose Q = ±1. The parameter e de- 
termines the proportionality coefficient between the po- 
tentials: Ve^t{r)/ g{r) = —e/Q. Then the change of the 
dependent variables ^(r, tlj{r]^(p^t)^ allows one to 
reduce Eq. ^ to the 2D form 

idt^ = g{7j, ip) (-^2 -dl-e + g\i>f) V-. (6) 

Here g{r]^ (p) is the positive definite function, defined as 
g{r]^(p) = |V7^p/2 = |V(^p/2, where the gradients must 
be expressed in terms of r] and (p (see the examples be- 
low). The obtained equation ([6]) is considered for > 
and has to be supplied with the zero boundary condition 
V^(7y = 0,(^,t) = 0. 

First, we concentrate on time-independent solutions 
of Eq. ([!]): ^(r, t) = tlj{r]^(p)^ where tlj{r]^(p) solves the 
2D stationary GP equation with constant coefficients 
si/j = —{d^ +9^)'0 + ^|?/^PV;. Particular solutions of this 
equation can be represented as 



where u is a. constant and the real-valued function (j){r]) 
solves the problem 

£4>{v) = -d^^4>{v) + 0\<Piv)f4>{v), m = o, (8) 

with E = e — v'^ and r] > 0. 

Turning to the physical meaning of the obtained so- 
lutions, we observe that it follows from the ansatz, 
ilj{r],(p) = e^^'^(j){r])^ that V(/:^(r) can be identified as the 
superfiuid velocity. Hence, the introduced analytic func- 
tion f{z) is nothing but the complex potential of the 
respective two-dimensional fiow. As it is well known [l2| 
such a potential defines the current, Jc of the fiuid 
through a given contour C C I), as well as the circu- 
lation Tc along C: 

[ f\z)dz = Vc^iJc (9) 
Jc 

(the prime stands for the derivative with respect to z). 
Since in our case f{z) is analytic, this integral is zero 
for any closed contour C bounding a simply connected 
domain. Thus the described fiow has neither sources nor 
vorticity in D. 

We also observe that if the change of variables im- 
plies growth of |V(/:^(r)| with r, then the physical mean- 
ing might have only solutions with densities decaying at 
the infinity (thus ensuring decaying currents). Whenever 
one concerns with finite densities at the infinity, the phys- 
ically meaningful solutions would correspond to z/ = 0. 
This last constraint is assumed in what follows. 

Finally, we notice that if the contour dD included in 
the proposed scheme is closed, this implies that linear and 
(or) nonlinear potentials are divergent at some point (s) 
of the boundary. This is clear from the nature of the 
conformal mapping, since in this case there should be a 
point of the boundary, which is mapped into the infinity 
point. In its turn such a point gives an origin to the 
singularity of the linear and nonlinear potentials. Such 
cases will be exclude in what follows, although they still 
may have physical relevance. 

III. EXAMPLES OF EXACT SOLUTIONS 

While large diversity of the domains can be considered, 
here we concentrate on the two simplest representative 
examples for the conformal mapping (j2j). 

A. The quadrant x > \y\ with the boundary x = d=^ 

The first domain considered below, is given by Di = 
{x > \y\} (i.e., Di is the quadrant of the (x, ?/)— plane) 
with the boundary dDi = {x = y > 0}\J{x = —y^ y < 
0}. Respectively, the conformal mapping is chosen as 
f{z) = z^, and thus 77(r) = x'^ — y'^ and (p{r) = 2xy. The 
linear and nonlinear potentials are now given by 



(7) 



V;,t(r) = -2e\r\\ g{r) = 2g\r\\ (10) 
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i.e., they are the hnear expulsive parabohc potential and 
the parabolic nonlinearity. Now g{r]^(p) = 2^7^^ -\- (p'^. 

B. The strip < y < tt with the boundary 2/ = 0, tt 

Another domain explored below is a strip D2 = {x e 
< y < 7t} with the boundary dD2 = {x e y = 
0} U {x G M, ^ = tt}. Now the function performing con- 
formal mapping to the upper half-plane is f{z) = and 
the new variables are determined as 77 (r) = smy and 
(f{r) = cosy. The linear and nonlinear potentials al- 
lowing for the exact solutions are now given by 

Kxt(r) = -2£e2^ ^(r) = 2ae'^ (11) 

and g{r],if) = \ {r]^ ^(p^). 

C. Exact solutions 

Turning to the exact solutions, below we consider only 
the case of repulsive interactions ^ = 1, as the most 
natural candidate to produce stable stationary flows. As 
the two simplest solutions of Eq. ([1]) (see e.g., Refs. [6, 
fllj ) we study, the "dark soliton" shape 

^ds{v) = ^ds{v.v) = V^tanh|^y|r?(r)j e^'-^W (12) 
and the nonlinear periodic modulation 

with £ = e — v'^ > and /c G (0, 1] being the modulus 
of the Jacobi elliptic sn-function [notice that '0(is(O, ip) = 
'0sn(O, = 0, in conformity with the imposed boundary 
conditions]. These exact solutions given by Eqs. (fT2)) and 
(p!3|) and the corresponding velocity fields in the domains 
Di and D2 are illustrated in Fig. [TJ 

The obtained exact solutions, however leave several 
open questions related to their practical feasibility. First, 
the stability of the flows was not investigated, so far. 
Second, in all considered cases infinitely growing poten- 
tials were used, while any cut-off (which exists in the 
real world) may strongly perturb, and even destroy the 
solutions. To address these issues we now turn to direct 
numerical simulations. 

In the case at hand the potentials 14xt(i*) and ^(r) 
grow with X, while the density goes asymptotically to 
(or is bounded by) a certain constant level. Bearing this 
in mind we construct a physical system with bounded 
potential in the following way. For x < we consider 
shifted linear and nonlinear potentials 14xt(r — 1*0) and 
^(r— ro), where fq = (xq, 0) with xq being a constant shift 
vector, while lext(r) and ^(r) are given by the analytical 
formulas ([5]). To define the physical potentials at x > 
we mirror them at x = 0, thus obtaining ]4xt(— r+ro) and 
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FIG. 1. (color online). The profiles for exact solutions given 
by Eqs. (|12p and (|13p and the corresponding velocity fields 
in the domains Di and D2. Left column: Panels (a), (b), 
(c) show the velocity field V(/?(r) = (2^,2x), and densities 
|^ds(r)|^ and |^sn(r)|^ for solutions given by ^ and {[l3| 
in the domain Di. Right column: Panels (d), (e), (f) show 
the velocity field V(/?(r) = (e^ cos 2/, sin^), and densities 
\^ds{Y)\^ and |^sn(r)|^ for solutions given by ()12p and ()13|) 
in the domain D2. The parameters are e = 2, zy = 0, and 
k = 0.8. 

^(— r + ro). To confine the condensate along ^/-coordinate 
we introduce additional linear tarp potential as follows: 
Vadd = for r G L> and Vadd = Vb, where Vq is large 
enough, for r ^ I^. For more details see the next chapter 
where we discuss the numerical studies of the problem. 



IV. DIRECT NUMERICAL SIMULATIONS 

Now we can perform direct numerical simulations of 
Eq. ([1]) with the additional confining potential Vadd tak- 
ing the initial distribution of the field in the form given by 
the analytical formulas ^^(r — ro) for x < and V^(ro — r) 
for X > 0. As it is clear our ansatz does not satisfy 
Eq. ([1]) only along the line x = and in the areas ^ < 
and y > 7T. Since, however we impose strong confining 
potential, we expect that the field is very weak outside 
the stripe < y < ir and our ansatz stays sufficiently 
close to the real stationary solution of Eq. ([T]) with the 
introduced physical potentials. The analytically found 
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FIG. 2. (color online). The ansatz produced from the solu- 
tions the with parameters £ = 2, zy = 0. In the upper part 
of the figure the distribution of the linear potential along y 
is shown for x — 0. On the right the distributions of the 
linear and nonlinear potentials are shown alongside with the 
introduced loss for ^ = 0. In the areas marked by "1" and 
"2" the anzatz exactly coincides with the analytical solution 
(neglecting exponentially weak loss). 



solution and the corresponding potentials are shown in 
Fig. O Within the areas "1" and "2" the ansatzes used 
in the numerics as the initial condition coincides with the 
analytical solutions. However on the boundary between 
these areas the ansatzes do not satisfy the equation for 
stationary fields. The ansatz does not satisfy the equa- 
tion outside the areas either. However the density of the 
condensate is low outside the areas and so there are some 
reasons to believe that the ansatz is close to the station- 
ary solution. Let us remind here that in our numerical 
simulations to keep the condensate localized within the 
stripe we used strong linear potential, see the distribution 
of the Hnear potential along y. 

We carried out numerical simulations for the solutions 
obtained in the domain D2 and found out that some 
non-stationary excitations appear but the solution sur- 
vive and stay rather close to the initial field distribu- 
tion. To check the stability of the solution for longer 
time we had to get rid of the propagating excitations. 
Since we perform numerical simulations in finite win- 
dows the only way to eliminate the propagating excita- 
tions is to introduce losses in the area where initially 
the density of the condensate is exponentially weak. 
To do this we introduce linear losses in the form 7 = 
70 [exp(-(x - xiY /wl) + exp(-(x + xiY)/wl] . Then 
the losses are negligible in the area where our analyti- 
cal solution is big and so the solution is practically un- 
affected by the artificial losses used in the numerics. In 
the same time all the perturbations propagating away are 
quickly absorbed. In this case the excitations disappear 



(a) (d) 




FIG. 3. (color online). Panels (a) and (d) show the initial 
distribution of the density of the condensates for cases 1 and 
2; panels (b) and (e) show the distributions of the condensate 
densities at t = 100; panels (c)and (f) show the differences 
between the initial density distributions and the density dis- 
tributions at t = 100. The shifts of the initial conditions and 
the potentials are xq = —2.5 for (a)-(c) and xo = —2 for (d)- 
(f). The solution parameters are s = 2, zy = 0. The field 
is kept within the stripe by strong repelling linear potential 
Va = 10^ for y < and y > tt. The parameters for the linear 
losses are jo = 5, wo = 1, xi — 10. 

very quickly and one can see in Fig. [3] that the surviv- 
ing solution is very close to the initial distribution. As 
it is expected the stationary solution deviates from the 
ansatz most strongly around the points {x = 0^ y = 0) 
and (x = 0, ^ = tt). 

In the case of attractive interactions (of a negative scat- 
tering length) the solution appears to be unstable against 
collapse. 



V. THE GENERALIZED CONFORMAL 
MAPPING 

Next we address possibilities of getting more general 
types of the potentials, allowing for exact solutions. To 
this end we introduce the generalized relations [c.f. the 
Cauchy-Riemann equations (j3j)]: 

p'^{fi)d^fi{r) = dy(p{r), {f})dyf}{Y) = -d^(p{r), (14) 

where ^^(t)) is a positive-definite function of 7}(r) only 
[as it is clear the arguments presented below can be also 
applied to the case where p'^{(f) is a function of (f{r) 
only]. These equations still define transformation to the 
orthogonal coordinates (7), (^), however now satisfying the 
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relations [c.f. Eq. (j4])] 



linear and nonlinear potentials given by [c.f. Eq. (pTj) ] 



V-[p^(7})V7)] =0, V-[p^(7))V(^] =0, V7}-V(^ = 0. (15) 

Moreover, now we have that |V(^(r)p = p^(7))|V7)(r)p. 
Notice that the generalized case [c.f. Eq. ([15])] can be 
reduced to the case mentioned above [c.f. Eq. (|4])] in the 
special case p^{f}) = 1. 

Next we repeat the steps described above for the con- 
formal mapping, and by the direct algebra show that the 
generalized ansatz 



^/>(7?,(^)=p(^)e-^W,^[7?(r)] 



(16) 



solves the 2D stationary GP equation ([T]) provided that 
(j){r) satisfy the stationary equation (j8]) with ri{r) substi- 
tuted by f)(r) and the hnear and nonhnear potentials are 
given by 



^^ext(r) 



5(r) 



(17) 



Notice that now the linear and nonlinear potentials are 
not proportional to each other any more [c.f. Eq. (j5j)]. 
To construct a particular example, we choose 



fi{r) = 2r]{r)- 



, (^(r) = (p(r), p{fi)-- 



(18) 



with 7^(r) and (p{r) solving Eq. (j4]). One can ensure that 
this choice also satisfies the conditions ([15]). Then the 
simplest solution of Eq. ([8]) with the repulsive nonlin- 
earity {Q = 1) is given by we study, the "dark soliton" 
shape ([T2]) with 7^(r) substituted by 7}(r) and valid for 
the positive chemical potential S = e — u'^ > 0. 

Then, according to the generalized ansatz with 7)(r), 
(^(r), and p{fj) given by Eq. (fT8|) . as well as 7^(r) and (p{r) 
defined in the domain we obtain the "dark soliton" 



^ds{r) =VSp{fi) tanh 



f){r)\ex.p{iiye^cosy)^ (19) 



with T] = 77(r) = e^sin?/, which solves Eq. ([T]) with the 



Kxt(r) = 



a2x 



{3-lQ{4:£ + v^)-2ri[v^{3 + 2ri) 



8(2 + r?)2 

+2£:(4+7?)(7?2+4r/+8)]} 



g(r) = -e'-i2 + 7ir. 

The presented solutions found in the domain D2 be- 
have in a way similar to the previous example, see Fig. |3l 
The numerical simulations prove that the found analyti- 
cal solutions can indeed be used as good approximation 
for the stable stationary solutions of Eq. ([T]) describing 
systems with physically relevant potentials. 

VI. CONCLUSION 

To conclude, we have shown that exact analytical solu- 
tions can be obtained in a large class of two-dimensional 
Gross-Pitaevskii equations with inhomogeneous linear 
and nonlinear potentials, defined on bounded domains. 
The method of constructing the models is based on the 
properly defined conformal mapping of the given domain 
into a complex half-plane. In the context of applica- 
tions to Bose-Einstein condensates, the obtained solu- 
tions having nontrivial phase depending on spatial coor- 
dinates can be interpreted as superfiuid fiows. In the case 
of negative scattering length (repulsive interactions) the 
background fiows, i.e., ones having no zeros in the open 
spatial domain, appear to be stable. The obtained re- 
sults generalize previous studies devoted to construction 
of the exact solutions, using the self-similar transforma- 
tion, to the two-dimensional models given on bounded 
domains. Moreover, the ideas presented in this paper 
are also able to apply in two-dimensional cubic-quintic 
models, two-dimensional multi-component models, etc., 
and to design linear and nonlinear potentials for control 
of Bose-Einstein condensates and nonlinear optical fibers 
in limited spatial domains. 
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